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ABSTRACT 

We present the results of a numerical simulation of the history and future develop- 
ment of the Pleiades. This study builds on our previous one that established statistically 
the present-day structure of this system. Our simulation begins just after molecular 
cloud gas has been expelled by the embedded stars. We then follow, using an N-body 
code, the stellar dynamical evolution of the cluster to the present and beyond. Our ini- 
tial state is that which evolves, over the 125 Myr age of the cluster, to a configuration 
most closely matching the current one. 

We find that the original cluster, newly stripped of gas, already had a virial radius 
of 4 pc. This configuration was larger than most observed, embedded clusters. Over 
time, the cluster expanded further and the central surface density fell by about a factor 
of two. We attribute both effects to the liberation of energy from tightening binaries 
of short period. Indeed, the original binary fraction was close to unity. The ancient 
Pleiades also had significant mass segregation, which persists in the cluster today. 

In the future, the central density of the Pleiades will continue to fall. For the first 
few hundred Myr, the cluster as a whole will expand because of dynamical heating 
by binaries. The expansion process is aided by mass loss through stellar evolution, 
which weakens the system's gravitational binding. At later times, the Galactic tidal 
field begins to heavily deplete the cluster mass. It is believed that most open clusters 
are eventually destroyed by close passage of a giant molecular cloud. Barring that 
eventuality, the density falloff will continue for as long as 1 Gyr, by which time most of 
the cluster mass will have been tidally stripped away by the Galactic field. 

Subject headings: open clusters and associations: general, individual (Pleiades) — stel- 
lar dynamics — stars: formation — binaries: general 



1. Introduction 



Despite recent advances in the field of star formation, the origin of open clusters remains a 
mystery. It is now generally accepted that all stars are born within groups. These groups are at first 
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heavily embedded within molecular clouds, their members obscured optically by copi ous interstellar 
dust. By the time the stars are revealed, only about 10 percent are in open clusters (jMiller &; Scalo 
193; I Adams fc Mve^ l200ll The remainder are in either T- or OB associations, both destined to 
disperse within a few Myr. In contrast, the stars within ope n clusters a re gravitationally bound to 
each other, and the group can survive intact for several Gyr (|Friellll995l ). How do molecular clouds 
spawn these relatively rare but stable configurations? 

One intriguing aspect of the mystery is that open clusters are intermediate in their properties 
between T- and OB associations. The for mer are relatively sparse in projected stell ar density, 
and contain up to about 100 members (e.g., iKenvon &: Hartmannlll995l ; iLuhmanl 120071 ) . The lat- 
ter, as exemplified by the nearb y Orion Nebula Cluster, begin with extraordi narily high densit y 
(jMcCaughrean fe Staufferl Il994l ) and contain well over a thousand members ([Hillenbrand! 119971 ). 
far more than the eponymous O and B stars. A published compilation of Galactic open clusters 
(|MermiHiodlll995l ) shows them to have from a few hundred to roughly a thousand stars, i.e., just 
in the middle range. Apparently, systems born with either too low or too high a population and 
density are fragile, while the relative minority falling in between can survive. 

There is already an ext ensive literature on y oung, bound clusters, both observational and 
th eoretical (fo r a rev iew, see lElmegreen et al.ll2000l ). Models for their origin, dating back at least 



to ILada et al.l (| 19841 ) , have focused on the need for a high star formation efficiency in the parent 
cloud. A standard computational technique, using N-body simulations, is to create stars in a 
background potential well, remove that potential through various pr escriptions, and then assess the 



result (e.g., iGoodwin &; Bastian 



2006 



Baumgardt fc Kroupal 120071 ) . Some researchers using this 



the bound remnants of expanding OB associations (jAdams 



2000 



approach, implemented either analytically or numeri cally, have hypothesized that open clusters are 

In recent 



Kroupa et al. 



2001 



years, most theor etical ideas have been motivated by fluid dynami c al simulations of turbulent 



collapsing clouds (jKlessen et al 



2000 



Vazquez-Semademi et al 



2003 



Li et all 12003 ). While much 



insight has been gained from these collective investigations, there has generally been too little 
contact of the theory with actual groups. A clear advance would be made if we could establish 
empirically the original state of one or more observed clusters. We would then be in a position to 
gauge how these particular systems were produced by star formation activity in their parent clouds. 



As a first step in this direction, IConverse &: Stahlerl (|2008l . hereafter Paper I) undertook a 
quantitative study of the present-day structure of a well-studied, relatively nearby open cluster, 
the Pleiades. We derived statistically, using a maximum likelihood analysis, such key properties as 
the stellar density distribution, mass function, overall binary fraction, and correlation between the 
component masses of these binaries. The point of that study was to provide the endpoint for any 
calculation of the system's previous evolution. 

We now take the second step. The age of the Pleiad es has been determined, from observations 
of lithium depletion, to be 12 5 Myr (jStauffer et al.lll998l ). Using the publicly available code Starlab 



(jPortegies-Zwart et al 



2001 



, Appendix B), we have run a suite of N-body calculations over just 
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this time period, to find that initial state which evolved to the current cluster, as gauged by our 
previous investigation. In doing so, we also establish the detailed history of the group over that 
epoch, and even into the future. 

A key assumption here is that the Pleiades divested itself of cloud gas relatively soon after its 
birth. There is currently no direct means to assess the duration of the initial, embedded phase, 
either in the Pleiades or any other open cluster We may take a clue from T associations, which 
are still surrounded (but not completely obscured) by molecul ar gas. No s ystems are observed 



with ages exceeding about 5 Myr, a striking fact first noted by iHerbiel (|1978l ). Presumably, older 
groups consisting of post-T Tauri stars have already driven away their clouds and are merged 
observationally into the field population. If a similar embedded period held for the Pleiades, it 
indeed represents a small fraction of the total age. Hence, we can establish, with some confidence, 
the cluster's structure just after cloud dispersal. A future study will investigate, using a combination 
of gaseous and stellar dynamics, how this early configuration itself arose. 

In Section 2 below, we describe in more detail our approach to the problem. We define the 
parameters characterizing both the initial configuration of the cluster and the evolved system. 
We then outline our strategy for finding the optimal initial state, i.e., the one whose descendent 
matches most closely the current Pleiades. Our actual numerical results are presented in Section 3. 
Here, we give the detailed properties of the inferred initial state. We also describe how the cluster 
changed up to the present, and how it will develop in the future. One of our key findings is that 
the cluster's evolution did not proceed in the classic manner associated with dynamical relaxation 



(jBinney &: Tremainelll987l . Chapter 8); we explore the origin of this discrepancy. Finally, Section 4 
discusses the implications of our findings on the earlier, embedded evolution of this, and other, 
open clusters. 



Method of Solution 



2.1. Initial Cluster Parameters 



2.1.1. Density and Velocity Distribution 
As in Paper I, we model the P leiades as a perfectly spheri cal system, although the cluster 



is observed to be slightly elongated (jRaboud &: M crmi iod 1998). This elongation seems to have 
been created by the tidal gravity of the Galaxy ( Wielen 19741 ). which would have exerted influence 
throughout the cluster's dynamical history. We assume that this modest tidal stretching had 
negligible effect on the internal evolution, and that relatively few stars were lost by tidal stripping 
over the Pleiades age. Thus we can safely ignore the associated Galactic potential. Similarly, we 
ignore mass loss through stellar evolution, which is negligible for our adopted mass function, over 
the 125 Myr age of the Pleiades. In Sections 3.2 and 3.3, we present simulations that include both 
effects 
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Returning to our standard runs, we further assume that the cluster was in virial equilibrium 
following expulsion of the gas. Any significant departure from equilibrium would be erased on a 
dynamical time scale, about 10 Myr for our input parameters and therefore much shorter than 



the evolutionary span of interest. Two popular choices for spherical equilibria are iKinel (|1966l ) 
models and polytropes. As will be explained in Section 2.3, King models do not include low enough 
density contrasts for a full exploration of initial states. We therefore used polytropes, which are 
more versatile in this regard. In polytropes, the stellar distribution function /, i.e., the number of 
stars per volume in configuration and velocity space, is given by 

J V ; [ £ < . V ' 

Here, A is a normalization constant, while £, the relative energy per unit mass, is 

£ = #(r) - v 2 /2 . (2) 

In this last equation, v denotes the stellar speed. Thus, £ is the negative of the physical energy, 
and also has an offset in its zero point, as conventionally defined. This offset is embedded in the 
relative potential \£(r), which is related to the usual gravitational potential <3?(r) by 

¥(r) = $(r t ) - $(r) . (3) 

The tidal radius r% marks the outer boundary reached by cluster stars. By construction, the relative 
potential Vl/ is positive inside the cluster and falls to zero at r = rt- 

The number of stars per unit volume is found by integrating the distributio n function over 



velocity space. The manipulations here are standard (IBinney fe Tremaind 119871 . Chapter 4), so 
we give only the essential results. We let m denote the stellar mass, assumed provisionally to be 
identical for all cluster members. Then p, the mass density of stars, is 

/""max 

p(r) = AirmA £ n ~ 3/2 v 2 dv . (4) 



o 



Here, v max (r) = \J2 W(r) is the maximum speed for a star at radius r. For such a star, £ = 0. 
The total physical energy per unit mass, $(r) + v 2 /2, is so the star can just reach rt- Using 

equation (2), equation (4) becomes 

; /•vW/ 2 \ ™-3/2 

p(r) = AirmA^-^ 2 J ( 1 — ^-J v 2 dv . (5) 

We define a new variable 9 = arcsin (v/V2^/), so that 

^r,n fir/2 

p(r) = 2 5 / 2 vrmA — / cos 2n ~ 2 d6 (6) 
n Jo 

= (2irf'*mA T §-W tt"(r) ■ (7) 
I\n + 1) 
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To calculate the relative potential VP(r), we use Poisson's equation. For our spherical system, 

where po and \&o are the central values of p[r) and \&(r), respectively. We define a dimensionless 
potential as tp = ^>/^>q, and a dimensionless radius as £ = r/ro, where the scale radius ro is 



y 4 vr G po 

Since p = p ip n , the new potential obeys the Lane-Emden equation: 

hWf)=- r ' <io> 

with boundary conditions ip(0) = 1 and ^'(0) = 0. The nondimensional tidal radius £j = r^/ro is 
the point where ip falls to zeroQ 

For any chosen polytropic index n, equation (10) can readily be solved numerically. Our task 
is to translate this nondimensional solution into a physical model of the initial cluster. Given n, 
the basic quantities characterizing the cluster are: iVtoti the total number of stars; r v , the virial 
radius; and m, the mean stellar massH The virial radius is defined by 

GM 2 , , 

r v = . 11 

2W V ' 

Here, M = N tot m is the total cluster mass, while W is the gravitational potential energy: 

1 f rt 

W = - 4:Trr 2 p$dr . (12) 

2 Jo 

In Appendix A, we show how to obtain the dimensional quantities ro, po, and \&o from our three 
input parameters and the solution 

Given the scale factors ro and po, we know the dimensional mass density p(r). We then 
populate space with stars according to a normalized distribution pi(r) such that p\{r)dr is the 
probability a star is located between r and r + dr. This probability is simply 

4?rr 2 p 

Pi(r) = -jr= • (13) 



x The value of £t is derived within each polytropic model. We stress that, despite the nomenclature, this "tidal" 
radius bears no relation to the truncation created by the Galactic potential. In Section 3.3, we describe simulations 
that include the external field. 

2 Unlike jV to t and r v , the mean mass m is not an independent parameter, but follows from our specified mass 
function and prescription for binaries; see §2.1.2 below. 
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The actual position vector of each star is then distributed isotropically within each radial shell. 



Finally, we require an analogous distribution for stellar speeds. At a given radius, the speed 
must be consistent with the prescribed energy distribution. Let p2(v\r)dv be the probability that 
the speed lies between v and v + dv given that its radius is r. Clearly, 



P2(v\r)dv x pi(r)dr 



f d 3 vd 3 7 



tot 



Replacing d s r by Airr 2 dr and d 3 v by Attv 2 dv, we have, after using equation (13), 

Air v 2 m f 



P2(v\r) 



P 



We take / = f{£) from equation (1) and use the definition of £ from equation (2), finding 



P2\v\r) 



2 T(n + 1) 3/2 
vr T{n -1/2) 



1 



v 

2~¥ 



n-3/2 



(14) 



(15) 



(16) 



Here, the relative potential is calculated at each r from ^ = ?/>(£)> where we recall that £ is 
the nondimensional radius. Given the stellar speed, i.e., the magnitude of the velocity vector, the 
direction of that vector is again distributed isotropically in space. 



2.1.2. Stellar Masses: Single and Binary 



Thus far, we have described a cluster that is composed of members with identical mass. In 
actual practice, we assign masses to the stars according to a realistic distribution. The parameters 
of this mass function for the initial cluster are among those we vary to obtain an optimal match 
between the evolved system and the present-day Pleiades. In the course of evolution, some stars 
will be given enough energy, through three-body interactions, to escape the cluster. The most 
massive ones die out over 125 Myr. It is therefore not obvious that the initial mass function is 
identical to that found today. 

We suppose that the distribution of stellar masses in the young Pleiades was similar in form 
to the initial mass function for the field population. In recent years, large-scale surveys of low- 
luminosity objects, combined with sp ectroscopy, have e stablished an accurate initial mass function 
down to the brown dw arf limit (e.g., ICovey et al.ll2008l ). The consensus is that the original power 
law of ISalpeterl ()1955l ) for masses above solar is joined at the lower end by a log normal function . 
This basic form appears to hold in diverse environments, including young clusters (|Chabrierll2005l ) . 



Let <f>{m) dm be the probability that a star's mass (in solar units) is between m and m + dm. 
We posit that this probability is 



(B/m) exp (— y 2 ) 
Cm a 



m min < m < fl 

fj, < m < m ma? 



(17) 
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where B, C, a, and the joining mass \i are all constants. We set m m i n = 0.08 and m max = 10 
(although we also tested a higher mass limit; see §3.2 below). The variable y is given by 

logm - logm f . 

V = y=- • (18) 

Here, tuq is the centroid of the lognormal function, and a m its width. For input parameters a, mo, 
and cr m , the constants B, C, and fi are determined by the normalization condition 

(m) dm = 1 , (19) 

and by requiring that $>(m) and its first derivative be continuous at m = \i. Analytic expressions 
may be found for the three constants, which we do not display here. 

Most stars are not single objects, but have binary companions. Indeed, the Pleiades today is 
especially rich in binaries (see §3.4 of Paper I). Such pairing must have been present in the initial 
cluster. We therefore view the positions and velocities established in the previous subsection as 
pertaining to iVtot stellar systems, rather than individual stars. Similarly, the symbol m used, 
e.g., in equation (13), actually denotes the average system mass, after accounting for binaries. 
We specify the global binary fraction as a parameter b, which gives the probability that a system 
actually consists of two stars. Conversely, a fraction 1 — b of the systems are indeed single stars. 
Their mass is distributed according to the probability (f>(m). 

Our analysis of the present Pleiades in Paper I showed that the masses of the component 
stars within binaries are correlated. Such a correlation must also have been present at early times. 
Accordingly, we include the effect in our initial state. Within the fraction b of systems that are 
binaries, we first independently assign masses to each component, using the probability distribution 
(j)(m). After identifying the primary mass m p and the secondary mass m s , we then alter the latter 
to m' s , where 

/ \ 7 



m' s = m s [^\ . (20) 
\m s J 

Here, 7 is an input parameter that measures the degree of mass correlation within binaries (see also 
eq. (42) of Paper I). If 7 = 0, the component masses are uncorrelated, while 7 = 1 corresponds to 
perfect matching. 

We give our binaries randomly inclined orbital planes, and a period and eccentricity dist ribution 



characteristic of p resent, solar- type bi naries both in the field (jDuquennoy Mayorlll99ll ) and in 



;it y aist ri 
■mI Il99lh 



the Pleiades itself (jBouvier et al.lll997l ). If P3(V) dV is the fraction of systems with periods between 



V and V + dV , then pz{V) is lognormal: 

p 3 (V) = ^=^exp (-z 2 ) , (21) 

where 

z = l ° gV - l ° gV ° . (22) 
V 2 cr-p 
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We set the centroid period to log Vq = 4.8 and the width to a-p = 2.3, where the period is measured 
in days. The eccentricity distribution has a thermal distribution: 



p 4 (e) = 2e, 



(23) 



as motivated both by observations (jDuquennoy fc Mayorlll99ll ) and theory (lHeggidll975l ). 



The initial cluster described thus far is homogenous, in the sense that any volume containing an 
appreciable number of systems has the same average system mass. However, there have long been 
claims of o bserved mass segregation in young clusters , i.e., an increase of av e rage stellar mass t oward 



the center (jSagar et al.lll988l ; iJones &; Stauffer 



1991 



Moitinko et al. 



1997 



Stolte et al.1 12009 1. The 



present-day Pleiades also exhibits this phenomenon to a striking degree (see §4.2 of Paper I). We 
want to see if this property developed on its own or was inherited from an earlier epoch. We 
accordingly include a quantitative prescription for mass segregation in our initial state. 

One system has a higher probability of being near the cluster center than another, in a time- 
averaged sense, if its relative energy £ is greater. Mass segregation ther efore manifests itself as a 
correlation between the system mass m and £. This fact was noted by iBaumgardt et al.l ([2008), 
who used it to implement a specific procedure for mass segregation. Here we have adopted a variant 
of their method that allows us to include the effect to a variable degree. We first assign £- and 
m-values to all member systems according to equation (1) and our prescription for binary masses. 
We then place the systems in two lists - the first ordered by increasing m, and the second by 
increasing £ . When we first construct these lists, the ranking of the system in the first is unrelated 
to its ranking in the second. This is the case of zero mass segregation. There would be perfect 
mass segregation if the two rankings were identical. 

Let us quantify the intermediate case. For a star of given mass, we find its index in the mass- 
ordered list. To assign an energy to that star, we choose the second (energy-ordered) index from 
a Gaussian distribution centered on the mass index. The width of this distribution, denoted as, 
can be infinite (no mass segregation) or zero (perfect segregation). More generally, we define a 
parameter (3, the degree of mass segregation, which varies between and 1. After some trial and 
error, we adopted the following prescription relating the width as to /3: 

1 



a £ 



N tot ln/3 . 



(24) 



The logarithmic dependence on f3 ensures that as has the desired behavior in the extreme limits. 
The proportionality with jVtot ensures that our algorithm gives the same degree of biasing in clusters 
of any population. 

In summary, (3 becomes another input parameter that we vary within the initial configuration. 
As we will see, having a non-zero f3 is critical to obtaining a proper match between the evolved 
cluster and the Pleiades today. Mass segregation was therefore present at a relatively early epoch. 

Since relatively massive stars preferentially reside near the center, our imposition of mass 
segregation alters the shape of the gravitational potential from that of a single-mass polytrope. 



- 9 - 



Relative to the total gravitational energy, the total kinetic energy has a value slightly below that for 
virial equilibrium. We rescaled all stellar velocities by a uniform factor to restore exact equilibrium. 
In practice, this factor was typically about 1.05. 

For convenient reference, Table 1 lists the full set of our input parameters for the starting 
state. Anticipating the results detailed in §3 below, the last column gives the numerical value 
of each parameter in the optimal configuration. We also list the associated uncertainties. The 
meaning of these uncertainties and how they were assessed, will also be discussed presently. 



2.2. Characterizing the Evolved Cluster 

After evolving a particular initial state for 125 Myr, we compare the outcome with the actual 
Pleiades. In making this comparison, it is important to "observe" the simulated cluster under the 
same conditions as the real one. Thus, we project the three-dimensional distribu tion of stars onto 



a two -dimensional plane, assumed to lie at the mean Pleiades distance of 133 pc (jSoderblom et al 



2005). The angular separation AO between each pair of stars is then determined. If A9 ICS denotes 



the telescope resolution, then any pair with A9 < A# res is taken t o be an unresolved p oint source. 



For the near-infrared catalog of the Pleiades analyzed in Paper I (jStauffer et al.l 120071 ). an appro- 
priate value of A6 Tes is 10 arcsec. Note that our unresolved sources include a small fraction (less 
than 0.5 percent) of triples and high-order systems, as well as a few unrelated pairs observed to be 
close in projection. We denote as N s the total number of point sources out to a radius from the 
cluster center of 12.3 pc, corresponding in angle to 5.3°. This was the radius enclosing the catalog 
of sources used in Paper 1. For each simulated evolution, we compare the final iV s -value with the 
observed Pleiades figure. 

The vast majority of stars observed in the Pleiades today are on the main sequence (see Fig. 1 
of Paper 1). The number of post-main-sequence objects, while relatively small, is sensitive to the 
shape of the stellar mass function. Hence, it is important that we reproduce, as closely as possible, 
the number inferred for the present-day cluster. At an age of 125 Myr, the main-sequence turnoff 
is about 4 Mq. If N± denotes the number of stars (singles or primary stars in binaries) whose mass 
exceeds 4 Mq, then this quantity, as calculated, may also be compared directly with that in the 
Pleiades. Similarly, we compare M tot , the total mass of all stars in the evolved cluster, with the 
Pleiades mass obtained through the statistical analysis of Paper 1. 

One striking result of Paper 1 was the prevalence of binaries. Specifically, we found that the 
near-infrared fluxes of the catalogued point sources demanded that the fraction 6 unres = 0.68 were 
unresolved binaries. (Any resolved binaries were listed in the catalog as separate sources.) For our 
assumed resolution limit A6 res , we could also assess 6 unres computationally for each evolutionary 
run. Again, this is the fraction of point sources representing two or more unresolved stars. Note 
that 6 U nres is less than the initially imposed binary fraction b, both because some pairs are wide 
enough to be resolved, and because others are torn apart in the course of evolution. 
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The distribution of stellar masses is, of course, another property that should be compared 
with the actual Pleiades. As just described, we set the form of the distribution within the initial 
configuration as a lognormal function with a power-law tail. The apportionment of masses within 
the evolved state could in principle differ, due to the escape of some stars from the cluster and 
the death of others with sufficiently high mass. Our procedure is first to find, within the output 
state, the normalized distribution of single stars. Included in this distribution are both isolated 
stars and the components of resolved binaries. We then peer within unresolved binaries and find 
the analogous distributions of primary mass m p and secondary mass m s . Finally, we record the 
distribution of the binary mass ratio, q = m s /m p . 

In Paper I, we statistically determined the stellar masses from the photometric data by as- 
suming that the single-star distribution was a pure lognormal, with centroid mo and width a m . 
For consistency, we characterize the evolved cluster in our simulations in a similar fashion. Now 
for a given single-star function and binary correlation parameter 7, the primary, secondary, and 
^-distributions are all uniquely determined. Appendix B outlines the mathematical derivation. The 
task is to vary 7, as well as tjiq and o~ m , for the presumed lognormal single-star distribution until 
this function, as well as the primary, secondary, and ^-distributions, best fit those we find directly 
in the numerical output. We then compare 7, mo, and a m to these same quantities derived in a 
similar way for the observed Pleiades. 

We next consider the projected density profile. We divide the cluster into radial bins that 
match those used in the analysis of the Pleiades. Th e resulting surface density of stellar systems is 



then fit to the empirical prescription of iKingl (jl962l ): 

E(R) = k ( - 1 = 1 \ . (25) 

Here, R is the projected radius, A; is a constant with the dimensions of a surface density, and R c 
and Rt are the core and tidal radii, respectively. We determine the values of k, R c and Rt which 
best match the data, i.e., the same parameters determined for the real Pleiades by an analogous 
fitting procedure. However, only R c is used in the final optimization routine (see below). We also 
determine the King concentration parameter ck = log (Rt/R c ), both for each evolved simulation 
and in the real cluster. From equation (25), the central surface density So is 

E Q = k I 1 1 ) . (26) 

^ Vi + (iW) v 

This is also compared to the Pleiades value. 

Finally, we measure the degree of mass segregation. For the evolved cluster, we compute the 
cumulative fraction of systems contained within a projected radius R, both by number (f^(R)) 
and mass (/m(-R))- As in Paper I (§4.2), the Gini coefficient is computed as 

G = 2 f\f M - f N ) df N . (27) 
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and then compared to that found in the Pleiades. 

Table 2 gives the full list of quantities evaluated for each evolved cluster. We also display the 
values found when the optimal initial state is used, as well as the corresponding figures in the actual 
Pleiades. Notice that mo and a m for the mass function do not match those in the initial state, 
as given in Table 1. This discrepancy arises partly from real changes of the stellar masses, but 
even more from our adoption of a simple lognormal when fitting the evolved cluster. The tabulated 
errors for the calculated quantities were obtained by running 25 simulations, all with identical input 
parameters. Since we populated the cluster stochastically, according to probability distributions 
(e.g., <j){m) in eq. (17)), initial states differed from one another in detail. The errors represent the 
standard deviations for each quantity in the evolved cluster, due solely to differing realizations of 
the initial state. The tabulated errors for the observed Pleiades are from the calculation of Paper I. 
In addition, N s and N4 are assumed to be Poisson-distributed, so that the errors are y/Ng and 
y/Nl, respectively. 



2.3. Optimization Procedure 

As a first guess, we set all the input parameters equal to values appropriate for the Pleiades 
today. In Paper I, we characterized the pr esent-day cl uster as a King model, so we initially adopted 



this prescription. The analytic models of iKind (|1966l ) have the distribution function 



f ^-{ £<0. (28) 

Here ^fo, the central value of the relative potential, is again set by our basic input quantities. The 
dimensionless parameter Wo, like ck, characterizes the degree of central concentration!! Running 
a King model for 125 Myr, we found that it invariably became more centrally concentrated than 
the actual Pleiades. We therefore tried successively lower Wo-values for the initial state. Now 
equation (28) shows that /(£) is proportional to £ in the limit of small Wo- Comparison with 
equation (1) reveals that such a model is equivalent to a polytrope of n = 5/2. Our search for 
low-concentration initial states therefore led us naturally to the polytropic models described in 
§2.1. Within the regime of polytropes with small n, we varied other input parameters, such as r v , 
as necessary. 

Once our computed cluster began to resemble the Pleiades, we changed to a more systematic 
gradient method for refining the initial state. Let x be the vector whose elements are the 9 input 
parameters listed in Table 1. Similarly, let y represent the 11 evolved cluster properties of Table 2. 
This latter vector is, of course, a function of x. To move y toward the values characterizing today's 
Pleiades, we need to evaluate, in some sense, the gradient of this function. 



3 The relation of Wo to ck is shown in Figure (4-10) of 
their plot, called \&o/o" 2 , is precisely Wo- 



Binnev fc Tremainel (1987). The independent variable in 
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A practical complication is one to which we alluded earlier. Even among evolutionary runs 
assuming an identical input vector x, the resulting y differs because of the stochastic sampling of 
the various assumed distribution functions. In computing the gradient, we need to take a step size 
h large enough that the resulting change in y exceeds that due to this re alization variance . We 



found that the prescription h = 0.5 x sufficed for this purpose (see §5.7 of lPress et al J 12002 . for a 
more rigorous justification). 

For each x, we first do 9 runs and average the result to obtain y(x). We then decrease, in turn, 
each element Xj to Xj — hj, and find the average output of two runs at each decreased x.,- value. 
Similarly, we find the average result of two runs at each Xj + hj. We thus establish the 11 x 9 
matrix of derivatives T>, whose elements are 

l iJ ~ 2h~ ■ ! - • 

The change in outputs for any subsequent input change Ax may then be approximated by 

Ay = VAx . (30) 

Here, the vector Ay is taken to be the difference between the current ^/-vector and that for the 
Pleiades. We may evaluate the 9 elements Axj by solving the 11 linear equations summarized 
in (29). Since the system is overdetermined, we did a least-squares fit to find that set of Axj which 
best satisfied the equations. 

As we took a step in x, we evaluated how close the resulting y was to y p , the aggregate 
properties of the observed Pleiades. We did a x 2 -test, where 

x2 = f : «y,>-y P , i )\ (31) 

1=1 a i 

Here, each < yi > is the average yi value, established by doing 9 runs with identical input values. The 
standard deviation <7j includes errors in both the inferred Pleiades properties and those generated 
by different statistical realizations of the input state: 

°i = <*l,i + °< yi > ■ ( 32 ) 



The first righthand term is the Pleiades variance whose square root is the error given in the last 

T <Vi> 



column of Table 2. The quantity <7< w > is the error in the mean y^. This error in the mean is related 



to a y i, the variance in each individual yi, by 



a %> = \ a li ■ ( 33 ) 



For the first few a:-steps, x 2 declined, but then stalled. Beyond this point, the gradient 
method itself was clearly failing, as it indicated initial states which evolved to configurations less 
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resembling the Pleiades. The difficulty was that the numerical derivatives of equation (29) were 
too crude to refine the initial state further. Refinements are possible in principle, but prohibitive 
computationally. After pushing the method to its limit, we were forced to stop the search before 
X 2 reached a true minimum. We took the last state in the sequence where x 2 declined to be the 
best-fit initial configuration. 

Our final task was to assess the errors in all input parameters for this state. These should 
reflect uncertainties in properties of the actual Pleiades, as well as the variation in output parameters 
among different runs using identical inputs. This latter effect is quantified by the covariance matrix 
y, whose elements are 

Vij = ((Vi- <Vi>) < Vj >)> • (34) 
The averagin g here refers t o different realizations using identical input parameters. Standard error 



propagation (jCowanl Il998l . Section 1.6) dictates that the known y is related to X, the desired 



covariance matrix of input parameters, through the derivative matrix and its transpose: 

y = VXV T . (35) 

We need to invert this equation to obtain X. As noted, the input errors should also reflect the 
observational uncertainties in the Pleiades itself. We do not know the correlation of these obser- 
vational uncertainties. Thus, we use on the lefthand side of equation (35) a matrix y' , formed by 
adding a 2 ^ to each diagonal element 3^- 

Since T> is not a square matrix, a standard inverse cannot be defined. Howe ver, the produc t 



is square, and so has an inverse, provided it is not singular. As discussed in iGravbilll (|1983l ). 
this fact allows us to define the pseudo-inverse of T>: 

V + = (V T Vy l V T . (36) 

The term "pseudo-inverse" is appropriate since 

V + V = (V T Vy l V T V = 1 , (37) 

where I is the identity matrix. Taking the transpose of this last equation, we also find 

(V + V) T = V T (V + ) T = 1 T = 1 . (38) 

By employing equations (37) and (38), inversion of the modified equation (35) is straightforward: 

v + y' (V + ) T = V + VXV T (v + ) T 

= X . (39) 



The errors in the initial cluster parameters of Table 1 are then the standard deviations obtained 
from the diagonal elements of X. 
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3. Numerical Results 

3.1. Global Properties of the Cluster 

Table 1 lists the optimal values for the parameters characterizing the initial state. The poly- 
tropic index n is about 3, corresponding to a volumetr ic, center-tq -average, number density contrast 
of 540 This particular polytrope closely resembles a iKingj (|1966l ) model with Wq ~ 1.4. Note the 



relatively large uncertainty in the optimal n, reflecting the fact that a range in initial density con- 
trasts relaxes to a similar state after 125 Myr. There is much less uncertainty in the virial radius 
r v , which is surprisingly large compared to observed embedded clusters (see §4 below). Smaller 
assumed revalues, however, evolved to systems with too high a density contrast. 

Figure 1 shows, as the dashed curve, the initial surface density as a function of projected 
radius. Also plotted (solid curve) is the evolved surface density, along with observed data from the 
Pleiades. Notice how the surface density decreases with time, a result of the inflation experienced by 
the entire cluste r. This behavior contrasts with expectations from the standard acount of dynamical 



relaxation (e.g. iBinney Tremaind 119871 . Chapter 8). The swelling of the central region that we 
find is consi stent, however, with previou s simulations of binary-rich clusters with relatively low 
populations ( Portegies-Zwart et al.1 [200ll ) . We explore further the underlying physical mechanism 



in §3.3 below. 

Note from Table 1 that iVtotj the initial number of stellar systems, is determined to within 
about 5 percent uncertainty. The main constraint here is the need to match N s , the final, observed 
number of point sources. Note also that iVtot < N s throughout the evolution. Almost all the stellar 
systems are binaries. Some of these are wide enough that they could be resolved observationally. 
Thus, the total number of point-like (i.e., unresolved) sources is always higher than N tot , the 
number of stellar systems (resolved or unresolved). By the same token, the unresolved binary 
fraction, b unres = 0.68, is significantly less than the full initial binary fraction, b = 0.95. Indeed, 
we were forced to choose a 6-value close to unity in order to make 6 un res close to the observationally 
inferred figure (see Table 2). 

Figure 2 quantifies the degree of mass segregation in the evolved cluster, Following the tech- 
nique introduced in Paper I, we plot /m, the fractional cumulative mass at any projected radius, 
against f^, the fractional cumulative number. The fact that this curve rises above the dashed di- 
agonal (/a/ = fx) indicates the existence of mass segregation. The empirical Jm — In relation for 
the Pleiades, shown by the points with error bars, is well matched by the simulation. We were able 
to obtain this match only by adopting a non-zero value of /3, the initial degree of mass segregation 
defined in equation (24) | 



Because of mass segregation, the center-to-average contrast in the volumetric mass density is higher, about 100. 

5 Table 1 lists, for convenience, a symmetrical error on the best-fit p. Although the lower bound is accurate, even 
higher values give acceptable results, due to the saturation of mass segregation described in §3.2 below. 
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Four quantities in Table 1, mo, cr m , a, and 7, concern the mass function. The number of 
stars escaping the cluster during its evolution is relatively small (on average, 280 of the 2400 
stars present initially). Because of this small loss, and because few members evolve off the main 
sequence, the initial and final mass functions are essentially identical, and all four parameters are 
highly constrained by the observations. Note, in particular, that the exponent a of the power-law 
tail directly influences N4, the observed number of massive stars. The binary correlation parameter 
7 is independent of the single-star mass function, but influences the primary, secondary, and q- 
distributions, as described in §2.2. Any substantial variation in 7 would alter the corresponding 
parameter obtained statistically for the observed cluster (see Fig. 5 of Paper I) . 

Figure 3 compares our evolved single-star mass function with the Pleiades. The solid curve is 
a lognormal fit to the simulation result, which is fully characterized by mo and a m in Table 2. The 
data points, along with error bars, represent the inferred single-star mass function for the Pleiades, 
obtained through the maximum likelihood analysis of Paper I. The agreement with the simulation 
is naturally poorest at the highest masses, since we modeled the output as a pure lognormal, in 
order to be consistent with the procedure adopted in Paper I. 

Finally, Figure 4 shows the initial distribution of the binary mass ratio q. We see how nearly 
equal-mass systems are strongly favored for our best-fit 7 of 0.73. In the simulations, this distri- 
bution evolves almost unchanged, and closely matches the one inferred for the Pleiades today. The 
figure also displays the g-distribution for the hypothetical case of 7 = 0. Such random pairing of 
stellar masses does not result in a flat curve, as one might expect. Instead, it reflects the character 
of the single-star mass function, which here is lognormal. As seen in Figure 4, the g-distribution 
for 7 = peaks at q = 0.34 and still vanishes as q approaches 0. 

3.2. Past Evolution 

We can now describe, based on our suite of simulations, the evolution of the cluster from its 
initial state to the present epoch. The main trend is an overall expansion of the system. This 
tendency is clear in Figure 5, which shows the variation in time of the virial radius, r v . After an 
initial drop, lasting about two crossing times (t CT oss = 10 Myr) the radius steadily swells, increasing 
by about 40 percent to the present. From the definition of r v in equation (11), we infer that the 
gravitational potential energy W is decreasing in absolute magnitude, i.e., the cluster is gradually 
becoming less bound. Note that we do not obtain r v by calculating W directly, but through fitting 
the cluster at each time to a King model, and then finding the appropriate r v for the best-fit model 
parameters. 

Figure 5 shows that R c , the projected core radius, displays similar behavior to r v . After the 
transient phase which again lasts about two crossing times, R c also swells, albeit more slowly. 
Analogous early adjustments are evident in other global quantities (see Figs. 6-8). This transient 
results from our implementation of mass segregation, which alters slightly the gravitational potential 
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(recall Section 2.1.2). Although the initial cluster is in virial equilibrium, the stellar distribution 
function is no longer a steady-state solution to the collsionless Boltzmann equation. Within the 
first two crossing times, the distribution readjusts to become such a solution. The core radius 
bounces, before settling to a value that subsequently evolves more gradually. 

Expansion of a cluster's ou ter halo is one mani f estati on of dynamical relaxation. However, 



application of equation (4-9) of iBinnev &: Tremaind (|1987l ). with In A = In (0.4 N), reveals that 



the relaxation time is 250 Myr, or about twice the Pleiades age. In addition, the inner cores of 
relaxing systems shrink, giving energy to the halo. The secular expansion of R c further indicates 
that we are not witnessing the usual effects of dynamical relaxation. Figure 6 provides yet another 
illustration of this point. Here, we see that the King concentration parameter ck remains virtually 
constant, again following an initial adjustment. Recall that ck = log (Rt/R c ), where Rt is the 
projected tidal radius. Thus, Rt and R c swell at about the same pace. 

The projected surface number density, S(J2), currently peaks strongly at R = (Fig. 1). The 
actually central value, however, previously declined from an even higher level. Figure 7 shows this 
gradual decline, which is consistent with the previously noted rise in R c . Thus, R c increases from 
1.6 to 2.2 pc over the period from t = 30 Myr to t = 125 Myr. Over the same interval, So falls by 
a factor of 0.50, which is close to (1.6/2.2) 2 . The number of systems in the core therefore remains 
virtually constant as the core itself expands. The volumetric number density similarly falls in the 
central region. 

In Paper I, we documented a strong degree of mass segregation in the current Pleiades, quan- 
tifying this property through the Gini coefficient. Another result of the current study is that G 
did not attain its current value through purely stellar dynamical evolution. As seen in the top 
curve of Figure 8, G(t) rose only slightly at first, and then remained nearly constant, even declining 
somewhat in the recent past. Initial states in which the parameter j3 was too low never attained 
the requisite degree of mass segregation. As an illustration, Figure 8 shows also the result from 
a single simulation using j3 = initially. The Gini coefficient does grow, but not by enough to 
match observations. We note, parenthetically, that G(t) exhibits oscillatory behavior over the a 
period that roughly matches the crossing time. These oscillations (unlike the initial readjustment) 
were washed out in the averaging procedure that produced the top curve in the figure. Finally, 
we remark that G(t) appears to saturate in time. We will return to this interesting phenomenon 
shortly. 

All the simulations we have described thus far ignored any effects of stellar evolution. We could 
afford this simplification because of the relatively small number of cluster members that would have 
evolved significantly over 125 Myr. However, the code Starlab does have the capability of tracking 
stellar evolution, including mass loss, through fitting formulae. As a check, we retained our usual 
maximum mass of 10 Mq and ran 25 simulations using the best-fit initial cluster parameters, but 
with stellar evolution included. The mass loss from relatively massive cluster members did not 
have a significant dynamical effect, and the endstate of the cluster was essentially identical. With 
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reference to Table 2, the only parameter that changed appreciably was N4, which fell. 



In more detail, the few stars above 7 Mq usually evolved to white dwarfs of approximately solar 
mass. On average, about 10 white dwarfs formed, of which 7 were the secondaries within binaries. 
Even the few that were primaries were faint relative to their main-sequence compa nions, and thus 



would be difficult to detect. Our findings are thus consistent with the observation of lFellhauer et al 



(|2003l ) that white dwarfs are generally rarer in open clusters than might be expected statistically 
from the initial mass function. 

Finally, we relaxed the upper mass limit in the single-star mass function and allowed the 
maximum mass to be arbitrarily large, according to the power law in equation (17). Choosing stars 
stochastically from this distribution yielded a few members with masses as high as 40 Mq. If we 
again allowed for stellar evolution and used our standard initial cluster parameters, the evolution 
did take a different turn. The very massive stars represented a significant fraction of the total cluster 
mass, and their death had a quantitative impact. As before, the cluster went through an initial 
adjustment, partially from the heavier stellar mass loss. The system then smoothly expanded, but 
at a faster pace. At 125 Myr, the projected core radius R c was 3.1 pc, or 1.5 times larger than that 
of the present-day Pleiades. Similarly, the central surface density So was a factor 0.58 lower. Had 
we begun with very massive stars in an initial state a factor of 1.5 smaller than our standard one, 
a closer match would have resulted. 

These results were instructive, if somewhat academic. In reality, stars more massive than 
about 10 Mq would have inflated HII regions so quickly as to ionize and disperse the parent 
molecular cloud forming the Pleiades. In order to retain even a remnant, gravitationally bound 
cl uster, the initial mem bership must have been very large, about 10,000 sta rs in the simu l ation s 
of iKroupa et al.1 (|200lh . We stress that even this figure is a lower bound, as iKroupa et all (|200ll ) 
assumed a star formation efficiency in th e parent cloud of 33 percent by mass. Such an efficiency 
is plausible within individual dense cores (jAlves. Lombardi. &: Ladall2007l ). but significantly higher 
than observational and the oretical estimates in cluster-forming clouds (e.g. iDuerr. Imhoff. &: Lada 



1982 



Huff fc Stabled 120071 ). 



Suppose we nevertheless adopt this scenario as a limiting case, and assume provisionally that 
the Pleiade s progenitor contained at least 10,000 individual stars. Such groups are rare. Equa- 
tion (39) of lMcKee Williams! (|1997l ) gives the birthrate of OB associations based on their popula- 
tion of supernova progenitors (m > 8). If we use our adopted initial mass function to estimate this 
population, then the birthrate of relevant OB associations is 0.09 Myr^ 1 kpc~ 2 . This is a factor 



between 5 and 8 sm aller than the total formation rate of open clusters (lAdams &: Myers 



2001 



Miller &: Scald Il978l ). It is unlikely, therefore, that formation through dispersing OB associations 
dominates, and we continue to use an upper mass limit of about 10 Mq for the Pleiades. 
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3.3. Future Evolution 

The same calculations that reconstruct the past history of the Pleiades may also be used to 
predict its de velopment far into the future. It is still believed, following the original proposal by 



Spitzerl ()1958l ). that most open clusters are eventually destroyed by the tidal gravitational field of 
passing interstellar clouds, now identified as giant molecular complexes. In this project, we do not 
attempt to model encounters with such external bodies. However, Starlab can follow the effects 
of the Galactic tidal field, both through imposition of the appropriate external potential and by 
adding a Coriolis force to individual systems. We switched on the Galactic field, in addition to 
stellar evolution, and followed the cluster from its initial state for a total of 1 G yr. While most 



open clusters do not survive this long, some do last up to several Gyr (|Friell ll995). Our simulation 
thus models at least a portion of the Pleiades' future evolution^ 

Up to the present cluster age of 125 Myr, adding the Galactic tidal field and stellar mass 
loss made very little difference in the evolution. Beyond this point, the cluster will continue the 
overall expansion that characterized it in the past. As seen in Figure 9, the central density S 
keeps declining. The falloff is roughly exponential, with an e-folding time of 400 Myr. This figure, 
along with Figures 10 and 11, show average results from the 4 runs we conducted. Even after 
averaging, the calculated S D displays increasing scatter for t > 700 Myr. By this time, the total 
population has also fallen to the point that numerical determination of S D (through a fitted King 
model) becomes problematic. 

The decline in the cluster population, which was modest until the present, accelerates as 
stars are tidally stripped by the Galactic gravitational field. It is the lighter stars that populate the 
cluster's outer halo and that preferentially escape. Consequently, the average mass of the remaining 
cluster members rises. Figure 10 shows both these trends. Displayed here is N s , the number of 
systems contained within the initial Jacobi radius of 14.3 pcQ By 1 Gyr, the total membership 
has fallen to a few dozen systems. Meanwhile, the average system mass (m) rises, almost doubling 
by the end. Careful inspection of Figure 10 shows that (to) initially fell slightly to its present-day 
value. This falloff reflects the loss, through stellar evolution, of the most massive members, an effect 
which is eventually overwhelmed by the escape of the lightest systems. Notice again the jitter in 
the (to) curve at later times. 

Despite this qualitative change in the cluster's internal constitution, the degree of mass segre- 
gation remains essentially constant until very late times. Figure 11 displays the Gini coefficient. In 
detail, G(t) exhibits oscillations qualitatively similar to what we saw in lower portion of Figure 8. 
Nevertheless, its average magnitude does not appreciable change until t ~ 700 Myr, when it be- 



6 The very oldest clusters have large galactocentric radii, and thus experience both a weaker tidal field and less 
frequent encounters with giant molecular clouds. Clearly, the Pleiades does not fall into this category. 

7 The Jacobi radius is t he spherical average of the zero- velocity surface in the presence of the Galactic tidal field 
( Binnev fc Tremaine|[l987l . p. 452). 
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gins a steep descent. The very large scatter during this late epoch again reflects the diminishing 
population. 

What accounts for these trends? During most of the evolution, some process is able to enforce 
mass segregation, despite the continual depletion of the lightest members. This process evidently 
loses its efficacy at late times, when the total population falls too low. Earlier, we showed that 
dynamical relaxation did not establish the present-day level of mass segregation. Although we are 
now spanning a period well in excess of the initial relaxation time (250 Myr), we still do not see the 
classic behavior - monotonic shrinking of the core that feeds halo expansion. Up to t ~ 200 Myr, 
the projected core radius R c continues the increase noted earlier. Between 200 and 600 Myr, when 
75 percent of cluster members escape, R c does decline slightly, from 2.2 to 1.2 pc. Thereafter, the 
core swells once more. 

We believe that the system's overall expansion is due principally to the release of energy in 
three-body encounters, specifically, close passages of binaries and single stars. Over the Gyr time 
span, mass loss through stellar evolution and tidal stripping weakens the cluster's gravitational 
binding, rendering it increasingly responsive to such internal heating. We ascribe the maintenance 
of mass segregation, i.e., the inward drift of more massive stars, to dynamical friction with the 
background population. We shall revisit these key processes, binary heating and dynamical friction, 
momentarily. 

Figure 12 shows graphically how the cluster will appear far in the future. Here we show 
positions of the member systems projected into the Galactic plane, both at the present time and 
at t = 700 Myr. One sees at present a slight elongation along the direction toward the Galacti c 



Center. This tidal stretching is well documented observationally (jRaboud fc Mermilliodl Il998l ). 
Stars that leave tend to do so along that direction. But any appreciable excursion leads, because 
of the Galaxy's differential rotation, to a change in angular speed. As a result, two tidal streams 
develop that are orthogonal to the Galactocentric radius. These streams are present in both panels 
of Figure 12, but are especially noticeable in the diminished cluster shown at the right. 

Since we suspected that binaries were important in the gross dynamics of the cluster, we 
recalculated the entire 1 Gyr evolution after effectively removing all primordial binaries from the 
system. We began with the same, best-fit initial state as before, which had the usual distribution 
of single-star masses and initial binary fraction b = 0.95. However, we replaced every binary by 
a single star, located at the system's center of mass and comprising the total of the component 
masses. In addition, we turned off both stellar evolution and the Galactic tidal field, in order to 
explore the evolution under the simplest conditions possible. Note that our procedure for fusing 
binaries into single stars preserved the total cluster mass, number of stellar systems, and mass 
distribution of those systems. In other words, all two-body interactions between cluster members 
were the same as before. The important difference is that we eliminated the source or sink of energy 
associated with the internal motion of binaries. 



The results were both surprising and illuminating. The cluster still undergoes overall expan- 
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sion. Figure 13 shows that the central density again falls steadily. The nominal e-folding time is 
again 400 Myr, but the actual decline is not well fit by an exponential. The root cause of the cluster 
expansion is that new binaries continually form and interact with other stars. This process occurs 
principally near the relatively dense cluster center, where the most massive stars reside, along with 
other, more representative members. The component masses in the new binaries are high, typically 
8 (m) . Prompt fo rmation of bin aries is a well-documented occurrence in systems initially containing 
only single st ars (lAarsethlll97ll ). and the formation rate is greatly enhanced at higher stellar mass 
(jHeggid Il975l ). In our simulations, only 3 or 4 of these systems exist at any time. Nevertheless, 
they are significant dynamically, because of the cluster's relatively low gravitational binding. 

Any such massive binary with a separation less than 8 x 10 4 AU = 0.4 pc is hard, i.e., has 
a gravitational potential energy exceeding the initial mean kinetic energy of all cluster members. 
Thus, even the relatively wide binaries formed in these simulations, with initial separa tions of order 
10 3 AU, a re capable of heating the cluster dynamically. As has long been appreciated (|Heggidll975l : 
Hutlll983l ). the encounter of a hard binary with a third star usually results in a harder (tightened) 
binary. Both this pair and the isolated star have more translational kinetic energy than before. 
The extra energy, which comes at the expense of the binary's tightening, is quickly transferred to 
other cluster members. 



The same dynamical heating operates, of course, in all stellar groups containing binaries. How- 
ever, very populous systems, such as globular clusters, have such high gravitational binding that 
almost all newly formed binaries are soft. In this case, energy exchange via three-body encoun- 
ters has a minor effect, and the classical picture of dynamical relaxation via two-body encounters 
applies. In relatively sparse systems like open clusters, both primordial and dynamically formed 
binaries inject so much energy that they impulsively change the velocity distribution function and 
qualitatively influence the course of ev olution. This st ochastic resetting of the velocities, which 
was emphasized in the classic study of iTerlevichl (j 19871 ). is a conspicuous feature of the Pleiades 
evolution, both past and future. 



4. Discussion 

While undertaken primarily to reconstruct the history of the Pleiades, our study has shed 
light on a well-documented, but still poorly understood, feature of open clusters generally - mass 
segregation. We demonstrated that the currrent, rather high degree of segregation in the Pleiades 
could not have been the result of dynamical relaxation from a pristine state with homogeneous 
mass distribution. First, the cluster has only been evolving for about half its initial relaxation 
time. Second, a hypothetical cluster starting with no mass segregation cannot reach the present 
level. Quantitatively, the Gini coefficient rises, but not enough (recall Figure 9). 

Two conclusions may be drawn. The ancient Pleiades must have already had substantial mass 
segregation before it drove off the gas. Some other process, unrelated to dynamical relaxation, 
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must drive this effect, and continues to do so long into the future (Figure 11). The most obvious 
candidate is dynamical friction. A relatively massive star moving through a lower-mass population, 
experiences a drag force, causing it sink toward the cluster center. The associated time scale fo r 
braking, tpy, can be substantially smaller than th e dynamical relaxation time i re iax (|Spitzerl ll969). 
According to iPortegies-Zwart McMillan! (|2002l ) , the quantitative relation is 



(m) 

3.3 ^relax 



(40) 



for a heavy star of mass m in a background of average mass (m). 



Portegies-Zwart &: McMillanl ([2002]) and other researchers have invoked dynamical friction to 



explain mass segregation, focusing on very popul ous clusters in which massive sta r infall leads to the 
runaway growth of a central black hole (see also lGiirken. Freitag. fc Rasioll2004l ). Our work reveals 
a further, curious aspect of the phenomenon. Figures 9 and 11 suggest, and further calculations 
confirm, that G(t) saturates, regardless of its initial value. Why does the degree of mass segregation 
level off? A possible explanation is that, as the most massive stars sink to the center, the population 
there becomes increasingly homogeneous. Since (m)/m rises, so does £df- I n other words, mass 
segregation through dynamical friction may be a self-limiting process. 

Returning to the prehistory of the Pleiades, another significant finding is the relatively large 
size of the initial state. The virial radius r v began at 4 pc, while the projected half-mass radius 
of the initial cluster was about 2 pc. For comparison, the observed half-light radii of embedded 
clust ers, as seen in the n ear infrared, range from about 0.5 to 1.0 pc, with some outliers on either 
side (ILada fc Ladall2003l ). Thus, the initial, gas-free Pleiades had a radius 2 to 4 times larger than 
typical embedded systems. It may plausibly be argued that the Pleiades is an especially populous 
open cluster, and therefore beg larger configuration, far outside the typical range. With 

this caveat in mind, our result suggests that the system expanded during its earliest, embedded 
phase. This swelling, which was accompanied, or even preceded, by mass segregation, could have 
been due to the loss of ambient gas during the formation process. Interesti ngly, observations of 
extra-Galactic clusters appear to show a similar, early expansion phase (see iBastian et al.l 120081 . 
and references therein). 

We have stressed the importance of binary heating to explain the global evolution of the 
Pleiades, both past and future. This is a three-body effect, not considered in classical studies of 
dynamical relaxation. As we indicated, binary heating is more effective in less populous systems, 
including open clusters. In the near future, we hope to explore further this general issue of stellar 
dynamics, i.e., the demarcation between systems that do and do not undergo classical, dynamical 
relaxation. This study will necessarily delve further into the role of binaries. We also intend to 
repeat our Pleiades analysis with another, relatively nearby system of comparable age, to ensure 
that the Pleiades results are representative for the entire class of open clusters. 



Steve McMillan, one of the authors of Starlab, provided crucial assistance throughout this 
project. Not only did he instruct us in the workings of the code, but he even debugged portions of 
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we thank James Graham and Chris McKee for their continued interest and provocative questions. 
S. S. was partially supported by NSF grant AST-0908573. 
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A. Physical Scale of the Polytropic Cluster 

As described in the text, the basic quantities characterizing the initial state are the polytropic 
index n, along with iVtot, r v an d m. We show here how to obtain from these the dimensional scale 
factors ro, po, and ^o- These scale factors, along with the dimensionless solution allow us to 

construct the physical cluster model, as described in the text. 

We first define a relative potential energy: 

1 f' rt 

W' = - 4irr 2 p^dr . (Al) 

2 Jo 

This has the same form as the true potential energy W, but uses the relative potential instead 
of the physical one <1>. Indeed, solving equation (3) for <& in terms of \& and subsituting into 
equation (12) for W yields 

W = -$(r t )M-W' (A2) 

- -^"W. (A3) 
2r t 

Since W is related to r v , we should next establish a relationship between W and other dimensional 
quantities. 



Following iKing] (jl.966), we define a nondimensional potential as 

>-2 jn+1 



P = / 4fff +1 d( (A4) 
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r'o Jo Po ^o 

2 W 



Po rl *o 

We also define a nondimensional cluster mass: 



(A5) 



H = I 47rC 2 ^ n ^ (A6) 
Jo 

_ 1 / A„J2 P 



., i Airr — dr 
r$ Jo Po 

M 
Po4 ' 

Both f3 and \x can be calculated using the solution Their ratio is 

P _ 2W 

2W'r fi 
4irGM 2 ' 



(A7) 

(A8) 
(A9) 
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where we have used equations (A7) and (9) to make the last transformation. Solving this equation 
for W' and substituting into equation (A3), we find 

GM 2 (\ 4vr/3\ IK s 

W = — l- + —ti) . A10 

2 \r t r /j, 2 J 

This last relation gives us more information about the virial radius. We now see that the 
nondimensional version, = r v /ro, obeys 



1 1 47T/3 

T = T + 



+ • (aid 



Thus, £ v can be obtained at once from the nondimensional solution. Since r v itself is an input, the 
dimensional scale radius, ro, can be obtained from 



Similarly, the central density po is 



while is found from 



ro = "f ■ (A12) 



M 

Po = ^3 (A13) 
^ . (A14) 



= 4ttGpo^ (A15) 
'47T&A GN tot m 
H ) r v 



B. Distribution of Component Masses within Binaries 



Let <f> be the normalized distribution of single star masses. In the text, we used the same 
functional notation when referring to the mass distribution for the initial cluster (see eq. (17)). We 
are now concerned with the evolved cluster. Our assumed functional form will be different, but we 
retain the notation for simplicity. We again let 7 be the binary mass correlation parameter, as in 
equation (20) of the text. Here we show how to find, for a given 7, the distribution of primary and 
secondary masses, as well as the distribution of the secondary-to-primary mass ratio q. Our final 
expressions for the various distributions are rather cumbersome and not especially illuminating; we 
therefore limit ourselves to outlining the derivation for a generic single-star function <p. 

We first let the primary and secondary masses have provisional masses m* and m* , respectively. 
Assume that both components within binaries are drawn independently from the same distribution 
4>. Then the two-dimensional mass function of the binaries is 



As explained in §2.2 of Paper I, the initial factor of 2 on the righthand side of equation (Bl) 
accounts for the different integration limits of $5 and <p (for the latter, see eq. (19)). 

To implement binary mass correlation, we consider new primary and secondary masses, m p 
and m s , related to the previous ones by 




Here, (m*, m*) Am* Am* is the probability of finding a system with primary mass between m* 
and m* + Am*, and secondary mass between m* and m* + Am*. This function is normalized so 
that 




m 



p 




(B3) 



m 



(B4) 



We are interested in the distribution function $b(m p ,m s ), which is 





(B5) 



After evaluating the Jacobian, we find 




) 



7/(1-7) 




) 



7/(1-7) 



$ fe (m p ,m s ) = 



1 - 7 \m p 



4>{m p ) (f) m s 



(B6) 



Let us first consider (j) p (m p ), the distribution of primary masses. This function is the integral 
of $(,(m p ,m s ) over all appropriate values of m s : 




(B7) 
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The largest mass a secondary can have, given the primary mass, is m p itself: 

m s , max = m p . (B8) 

However, the smallest mass is not m m i n . This is indeed the smallest mass for m*. The correlation 
of primary and secondary masses implies that the minimum for m s is 

— — ) • ( B9 ) 

Thus, the two integration limits in equation (B7) are themselves functions of m p . 

The secondary mass function, (f) s (m s ), is similarly found by integrating <f>b(m p ,m s ) over all 
possible primary masses: 

4> s (m s ) = I ® b (m p ,m s )dm p . (BIO) 

The smallest value a primary mass can be, for a given secondary, is m s : 

"i-p.min = rn s . (Bll) 

Somewhat surprisingly, the largest value is not necessarily m max , again because of the imposed 
correlation. To find the correct maximum, we solve equation (B4) for m*: 

m * = nlp = m J/7 ( m *)(7^)/7 . ( B12 ) 

Since 7 lies between and 1, the exponent of m* is negative. Thus, for a given m s , m p is greatest 
when m* is smallest. Since the lowest value of m* is m m i n , we have 

m p , max = m\^{m n ^- 1)h (B13) 
/ m \ 1/7 

= m min (— • (B14) 

V m min/ 

However, for m s > m m i n (m max /m m i n ) 7 , this equation says that m P!max > m min , which is impossi- 
ble. In summary, m Pjmax is given by 

(m s /m max ) 11 m s < m min 

"T-p.max — S , , s 7 l mD J 

J ^max ^ '"■min v^T-max/^minJ 

Finally, we need 4> q (q), the distribution of the binary mass ratio g = m s /m p . As a first step, 
we find the two-dimensional mass function &i,(m p ,q). Proceeding as before, we have 



®b{m p ,q) = $ b (m*,m*) 



d(m p ,m s ) 



(B16) 



3 (m p , g) 

= m p $ b (m p ,m s ) (B17) 

(B18) 



1-7 
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The desired distribution is the integral of 3>{,(m p , q) over suitable m p -values: 

/ 1 '^p,max(q) 

4> q (q) = / <f> b (mp,q)dm p . (B19) 

"' m p,min(q) 



As the notation indicates, the limits of m p are subject to the restriction of a fixed q. Now q itself 
is given in terms of m* (= m p ) and m* by 



bp 



(B20) 



7-1 



^ ) . (B21) 

\ m *sJ 

Solving the last equation for m p gives 

m p = m * s g -V(i-7) . (B22) 

Since the exponent of q is negative, and since q itself lies between and 1, we see that m p > m*. 
Thus, for any q- value, there is always some m* s for which Trip — wi max . We therefore set 

»™p,max(q) = m max • (B23) 

The smallest value of m p corresponds to m* = m mm . It follows that 

™ P ,min(q) = ™min . (B24) 



We may, in principle, perform the integrals in equations (B7), (BIO), and (B19) for any specified 
single-star function 4>(m). In practice, we choose a lognormal: 

<t>(m) = - exp (-y 2 ) , (B25) 
m 

where D is the normalization constant, and the variable y is given by equation (18) in the text. 
With this form of 4>(m), the integrations may all be done analytically, although we do not reproduce 
the rather lengthy results here. 
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Fig. 1. — Surface number density as a function of projected radius. The dashed curve represents 
our initial configuration, an n = 3 polytrope. The solid curve is a iKind (|1962| ) model fit to our 
simulation results for the evolved cluster. Table 2 lists the parameters for this optimal model. The 
numerical results displayed are an average of 25 simulation runs. Also shown are Pleiades data 
with error bars, taken from Paper I. 
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Fig. 2. — Fractional mass versus fractional number for the Pleiades. The solid curve shows the 
average results of our simulations. The crosses represent Pleiades data with error bars, taken from 
Paper I. The dashed diagonal line is the hypothetical result for zero mass segregation. 
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Fig. 3. — Single-star mass function for the evolved cluster. The solid curve is a lognormal fit to 
simulation data. Also shown are Pleiades data, with error bars, from Paper I. 
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Mass Ratio q 

Fig. 4. — Initial distribution of the mass ratio within binaries, q = m s /m p . The solid curve was 
obtained using a lognormal fit to the calculated single-star mass function, including the proper 
binary mass correlation parameter 7. The dashed curve is the hypothetical distribution obtained 
with the same single-star mass function, but with no mass correlation (7 = 0). 
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Fig. 5. — Evolution of characteristic radii. The upper curve shows the three-dimensional virial 
radius r v , and is an average over simulation runs. The lower curve shows the projected core radius 
R c , and is also an average. The data point in the lower right is the observed Pleiades value for R c , 
along with error bars. 
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Fig. 6. — Evolution of the King concentration parameter. The curve is an average over simulation 
runs. The data point is the observed Pleiades value. 



-37- 




Fig. 7. — Evolution of the central surface number density. Shown is the average over simulation 
results. The observed Pleiades value is represented by the data point. 
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Fig. 8. — Evolution of the Gini coefficient. The upper curve is an average over simulation results. 
To the right of this curve is the observed Pleiades value. The lower curve shows the result from a 
single simulation run in which the mass segregation parameter f3 was artificially set to zero. 
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Time t (Myr) 

Fig. 9. — Evolution of the central surface number density over a total time of 1 Gyr. The data 
point is the present-day Pleiades value, with errors indicated. 
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Time t (Myr) 

Fig. 10. — Evolution of the total number of stellar systems, N s , and the average system mass (m), 
over 1 Gyr. Both quantities refer to systems within the initial Jacobi radius of 14.4 pc. The data 
points show the current Pleiades values, with error bars. 




Fig. 11. — Evolution of the Gini coefficient over 1 Gyr. Note the large scatter at late times, 
reflecting the falloff in total cluster population. The data point to the left is the current Pleiades 
value, along with error bars. 
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Fig. 12. — Positions of Pleiades members projected onto the Galactic plane. The data are from 
a single, representative simulation, for the two epochs indicated. A terrestrial observer is located 
133 pc in the negative x-direction. The Galactic Center is in the same direction, but 8 kpc distant. 
Galactic rotation is in the positive y-direction. 
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Fig. 13. — Long-term evolution of the central number surface density for a cluster with no primordial 
binaries. The curve was obtained by averaging 9 simulation runs. The data point shows the current 
Pleiades central density, with errors. 
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Table 1: 


Initial Cluster Parameters 




Symbol 


Definition 


Optimal Value 


n 


polytropic index 


3.0 ± 1.3 


Ntot 


number of stellar systems 


1215 ± 59 


r v 


virial radius 


4.0 ±0.9 pc 


m 


centroid of mass function 


0.12 ±0.04 M Q 


Cm 


width of mass function 


0.33 ± 0.06 


a 


exponent in mass function 


-2.20 ±0.04 


b 


fraction of binaries 


0.95 ± 0.08 


1 


mass correlation in binaries 


0.73 ± 0.09 


P 


degree of mass segregation 


0.5 ±0.3 
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Table 2: Evolved Cluster Properties 



Symbol 


Definition 


Calculated Value 


Pleiades Value 


N s 


number of point sources 


1244 ± 32 


1256 ± 35 


N 4 


number of systems with m > 4 


13 ±4 


11 ±3 


Mtot 


cluster mass 


939 ± 30 M 


870 ± 35 M 


^unres 


unresolved binary fraction 


0.68 ± 0.02 


0.68 ± 0.02 


m 


centroid of mass function 


0.12 ±0.03 M Q 


0.14 ±0.05 M 




width of mass function 


0.49 ± 0.05 


0.46 ± 0.04 


1 


binary correlation index 


0.66 ±0.01 


0.65 ±0.05 


R c 


core radius 


2.2 ±0.4 pc 


2.0 ±0.1 pc 


c K 


King concentration parameter 


0.98 ± 0.09 


0.99 ± 0.04 


So 


central surface density 


36 ± 8 pc" 2 


40 ± 3 pc" 2 


G 


Gini coefficient 


0.18 ±0.02 


0.20 ± 0.02 



